%function [aa bb] = imageStackMedianLogGradients(M)
%  Given an image stack, x-size x y-size x 3 x N
%  return the median of that stack

function [medX medY] = imageStackMedianGradients(M)
if size(M,4) == 1
    M = permute(M,[1 2 4 3]);
end
sz = size(M);

% get relevant sizes
m1 = M(:,:,:,1);
N = sz(4);

%% iterative median
[medX, medY, countX, countY] = deal(zeros(size(m1)));

fprintf('Computing iterative gradient:+++++++++++!\n');
fprintf('Computing iterative gradient:');
for px = 7:-1:0
    countX = 0.* countX;  % reset the counters for this loop
    countY = 0.* countY;
    for ix = 1:N
        im = double(M(:,:,:,ix));   % get the gradient log image for next set
        [d1 d2] = gradient(im);
        
        bigX = find(d1(:)>medX(:));  % count if this is bigger or smaller than current gradient image
        countX(bigX) = countX(bigX) + 1;
        bigY = find(d2(:)>medY(:));
        countY(bigY) = countY(bigY) + 1;
    end
   
    % adjust X- median based on number of times pixels are more than
    % current median
    bigX = find(countX > ( N/2 + 1));
    littleX = find(countX < (N/2 - 1));
    medX(bigX) = medX(bigX) + 2.^px;
    medX(littleX) = medX(littleX) - 2.^px;
    
    % adjust Y- median based on number of times pixels are more than
    % current median
    bigY = find(countY > (N/2 + 1));
    littleY = find(countY < (N/2 - 1));
    medY(bigY) = medY(bigY) + 2.^px;
    medY(littleY) = medY(littleY) - 2.^px;
    fprintf('.');
end
fprintf('\n');